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ABSTRACT 

Accurate modeling of nonlinearities in the galaxy bispectrum, the Fourier transform of 
the galaxy three-point correlation function, is essential to fully exploit it as a cosmolog¬ 
ical probe. In this paper, we present numerical and theoretical challenges in modeling 
the nonlinear bispectrum. First, we test the robustness of the matter bispectrum mea¬ 
sured from A/'-body simulations using different initial conditions generators. We run 
a suite of Wbody simulations using the Zel’dovich approximation and second-order 
Lagrangian perturbation theory (2LPT) at different starting redshifts, and find that 
transients from initial decaying modes systematically reduce the nonlinearities in the 
matter bispectrum. To achieve 1% accuracy in the matter bispectrum at 2 : < 3 on 
scales k <1 /i/Mpc, 2LPT initial conditions generator with initial redshift 2 ; > 100 is 
required. We then compare various analytical formulas and empirical fitting functions 
for modeling the nonlinear matter bispectrum, and discuss the regimes for which each 
is valid. We find that the next-to-leading order (one-loop) correction from standard 
perturbation theory matches with Wbody results on quasi-linear scales for 2 ; > 1. 
We find that the fitting formula in Gil-Marm et al. (2012) accurately predicts the 
matter bispectrum for 2 ; < 1 on a wide range of scales, but at higher redshifts, the 
fitting formula given in Scoccimarro & Couchman (2001) gives the best agreement 
with measurements from Wbody simulations. 
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1 INTRODUCTION 

The large-scale structure of matter and galaxies in the Uni¬ 
verse is an excellent cosmological probe. Thus far, the two- 
point correlation function of galaxies has been widely em¬ 
ployed to study the large-scale structure, and broaden our 
understanding of the Universe by providing tight constraints 
on, e.g., the curvature of the Universe, initial conditions, and 
properties of dark energy (e.g. (Tegmark et al. 2006; Blake 
et al. 2011; Parkinson et al. 2012; Reid et al. 2012; de la Torre 
et al. 2013)). A complete description of non-Gaussian large- 
scale structure, however, must include higher-order statis¬ 
tics. Although the matter density fluctuations follow nearly 
Gaussian statistics at early times, the nonlinear process of 
gravitational instability generically drives the matter den¬ 
sity non-Gaussian. In addition, a nonlinear relation (bias) 
between the galaxy density and matter-density contrast and 
the distortion of the measured redshift of galaxies due to 
peculiar velocity further enhance the non-Gaussianity. The 


leading order statistic sensitive to such non-Gaussianity is 
the three-point correlation function, or its Fourier transform 
the bispectrum. 

The galaxy three-point function and bispectrum have 
been measured previously in IRAS (Scoccimarro et al. 2001), 
2dFGRS (Verde et al. 2002; Gaztahaga et al. 2005), SDSS 
(Nichol et al. 2006; Nishimichi et al. 2007; Marm 2011; 
McBride et al. 2011a,b; Gil-Marm et al. 2015), and Wig- 
gleZ (Marin et al. 2013), but have not been as fruitful as 
their two-point counterparts because theoretical modeling 
of higher-point statistics has not yet reached the level that 
is required for sophisticated data analysis. Future surveys 
promise orders of magnitude improvement in both galaxy 
number density and volume coverage, permitting better de¬ 
termination of the higher-order statistics of the cosmic den¬ 
sity field, and the huge resources that the community has 
put into those large surveys compels us to exploit the higher 
order statistics. Upon accessing all the cosmological infor- 


© 0000 The Authors 


2 N. McCullagh et al. 


mat ion, the three-point statistics will be a powerful probe 
of inflation (Jeong & Komatsu 2009; Sefusatti & Komatsu 
2007; Baldauf et ah 2011), nonlinear structure formation 
(Scoccimarro et ah 1998, 1999), and astrophysics such as 
galaxy formation (Baldauf et ah 2011). 

In order to extract information from the bispectrum of 
galaxies we must have an accurate model of the bispectrum 
of the underlying matter field. Cosmological A/’-body sim¬ 
ulations are the most direct way of studying the nonlinear 
evolution of the matter bispectrum. One consideration when 
modeling the bispectrum using simulations (in particular at 
high redshift) is transients from initial conditions (Scocci¬ 
marro 1998; Crocce et al. 2006). The initial positions and 
velocities of the iV-body particles are often generated by us¬ 
ing Lagrangian Perturbation Theory (LPT), either to first 
order (Zel’dovich approximation) or second order (2LPT). 
Because the initial conditions from the Zel’dovich approxi¬ 
mation and 2LPT include spurious decaying modes—which 
do not exist in the real Universe—on top of the desired 
fastest growing modes, one must wait a long enough time 
for the decaying modes to be sufficiently suppressed (Scoc¬ 
cimarro 1998). For the case of the power spectrum, where 
the Zel’dovich transients appear at non-linear order, 2LPT 
initial conditions with starting redshift of zini = 49 may 
be enough for capturing the nonlinearities at 2 : = 0 (Crocce 
et al. 2006). As we shall show below, however, correct simula¬ 
tion of the matter bispectrum at redshifts 1 < z < 6 requires 
even earlier starting redshifts with 2LPT initial conditions. 
This is because in the Zel’dovich approximation, the decay¬ 
ing modes in the matter bispectrum affect the leading order, 
and suppress the resulting matter bispectrum even on the 
very largest scales in the simulations. 

Parallel to numerical simulations is theoretical study of 
the nonlinear evolution of the matter density field. Standard 
perturbation theory (PT) (for a review: Bernardeau et al. 
(2002)) models nonlinearities in the matter density field as 
a pressureless, single fluid evolving under gravitational inter¬ 
action. In PT, the leading order, tree-level, matter bispec¬ 
trum can be calculated from the second order solution for the 
density field, giving us an expression that is valid on large 
scales. Beyond tree-level, one can include successive higher 
order corrections to the bispectrum, which may correctly 
model the non-linear bispectrum in the quasi-linear regime. 
Other approaches, such as resummed Lagrangian Perturba¬ 
tion Theory (LPT) have also been developed to model the 
bispectrum on quasi-linear scales (Rampf & Wong 2012). For 
even smaller scales, where the r.m.s. of the density contrast 
is of order unity, PT breaks down and the matter density 
field is in the fully nonlinear regime. Various phenomenolog¬ 
ical fitting formulas based on V-body simulations have been 
developed that claim to better model the nonlinear behav¬ 
ior of the bispectrum at low redshifts. There are two ver¬ 
sions of fitting formulas available in the literature: the first 
proposed by Scoccimarro V Couchman (2001) and more re¬ 
cently one proposed by Gil-Marm et al. (2012). In light of 
the recent development of renormalized perturbation theory 
(Crocce V Scoccimarro 2006), one can also try simplistic 
models of ‘renormalizing’ the power spectrum while keeping 
the tree-level vertex (the second order kernel) intact. That 
is, we assume that the nonlinear behavior of the bispectrum 
is completely described by the nonlinear power spectrum, 
so the usual tree-level expression is altered by replacing the 


linear power spectrum with the nonlinear power spectrum. 
To verify these formulas, we compare the three nonlinear 
models of the matter bispectrum with the results from N- 
body simulations with the smallest transient effect, and find 
the region of validity for each of model, in particular at high 
redshifts (1 < z < 6). 

In this paper, we shall address the aforementioned nu¬ 
merical and theoretical challenges on accurately modeling 
the matter bispectrum in the following order. In Sec. 2, 
we discuss various ways of visualizing the bispectrum, and 
introduce a flattening plot, which shows all of the mea¬ 
sured bispectrum on a range of scales, that we shall use 
throughout this paper. Sec. 3 reviews the standard pertur¬ 
bation theory approach to calculating the tree-level bispec¬ 
trum and leading-order transients from initial conditions. 
Sec. 4 presents the effect of transients on the measured bis¬ 
pectrum from a suite of iV-body simulations with different 
initial redshifts and initial condition generators (Zeldovich 
and 2LPT). We then discuss in Sec. 5 theoretical and empir¬ 
ical models of the nonlinear bispectrum. We show that the 
fitting formulae in Scoccimarro V Couchman (2001); Gil- 
Marm et al. (2012) accurately predict the bispectrum in dif¬ 
ferent regimes: Scoccimarro V Couchman (2001) is accurate 
on the widest range of scales and redshifts, and works partic¬ 
ularly well at high redshift (z > 1), whereas Gil-Marm et al. 
(2012) is the most accurate model at lower redshifts (z < 1), 
but over-predicts the bispectrum at high redshifts. We also 
compare these models to perturbation theory predictions for 
the bispectrum. We conclude in Sec. 6. 


2 VISUALIZATION 

The bispectrum B(A;i,fe,fe) of the Fourier space density 
contrast field S{k) is defined as: 

{S{ki)5ik2)S{k3)) = (27r)^B(fci,fc2,fc3)fe(fei + fe + ks) , 

( 1 ) 

where the is a Dirac-delta function. The bispectrum is 
non-zero only where the three wavevectors form a triangle 
due to statistical homogeneity of the matter density field. 
Furthermore, statistical isotropy dictates that the bispec¬ 
trum is independent of the orientation of the triangle, so 
the magnitudes of the wavevectors specify the bispectrum^. 
Without loss of generality, we impose the condition that 
ki > k 2 > ks throughout the paper. For later use, it is use¬ 
ful to consider various limiting triangular configurations for 
reference: squeezed {ki ~ /c 2 ^ fe), elongated (ki = /c 2 +fe), 
folded {ki = 2 /c 2 = Sfe), isosceles (/c 2 = fe), and equilateral 
(ki = fo = fe). For graphical illustration of each triangu¬ 
lar configuration, we refer the readers to Fig. 1 of (Jeong V 
Komatsu 2009). 

How should we visualize the bispectrum? Fig. 1 shows 
the leading order matter bispectrum as a three-dimensional 
color plot. As we shall show in the next section, the leading 
order matter bispectrum in PT is given by 

B{ki, k 2 , ks) = 2F^‘’\ki,k2)PL{ki)PL{k2) + (2 cyclic) (2) 

^ With redshift-space distortions, statistical istotropy is violated, 
and the bispectrum depends on the angle between the wavevectors 
and the line-of-sight direction. In this paper, we only discuss the 
matter density field in real space. 
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Figure 1. 3-D plot of the full tree-level matter bispectrum, 
with ki > k 2 > ks.The color reflects the magnitude of 
\og{B{ki,k2,k3)). 


where 


5 fei • fe / 1 1 \ 2 (fci • 

7^ 2 \kl^kl)^7 klkl 


( 3 ) 


is the second order PT kernel and PL{k) is the linear mat¬ 
ter power spectrum. The three axes of Fig. 1 represent the 
magnitudes of the three wave numbers, ki, k 2 , and ks in 
which the bispectrum is shown only for the tetrahedronic 
region defined by ki > k 2 > k^ and ki < k 2 -\- k^ (trianglu- 
lar condition). The amplitude of the bispectrum is coded by 
colors. We can see from Fig 1 that the bispectrum has the 
largest magnitude at around /ci ~ ~ 0.02 h/Mpc 

where the linear power spectrum peaks. We can also ex¬ 
amine the behavior of the bispectrum at the limiting trian¬ 
gular configurations, which are along the planes and lines 
that make up the shape of the region in which the bispec¬ 
trum is defined. This way of visualization is often used in 
the cosmic microwave background (CMB) literature when 
discussing the bispectrum of temperature anisotropies (Fer- 
gusson et al. 2012; Planck Collaboration et al. 2013). While 
the three-dimensional color plots are useful for showing the 
overall structure of the bispectrum, it is not apt for showing 
its detailed shape and scale dependence. 

A useful way of visualizing the shape and scale depen¬ 
dence of bispectrum is by plotting slices of constant ki . The 
plots in the left column of Fig 2 show slices of the bispec¬ 
trum in Fig. 1 at /ci = 0.01 /i/Mpc, ki = 0.05 /i/Mpc, and 
ki = 0.15 h/Mpc (from top to bottom). Note that the val¬ 
ues shown have been normalized by the maximum value of 
B(ki, k 2 , ks) in each slice. We also show the locations of the 
different triangular configurations in the same plot. From 
these plots, we can see that the squeezed configuration al¬ 
ways has a suppressed signal compared to other elongated 
triangles. On large scales, the bispectrum of equilateral tri¬ 
angles is higher than in other configurations, whereas on 
smaller scales, it is suppressed compared to, say, isosceles 
triangles. Jeong & Komatsu (2009) present a detailed dis¬ 
cussion about the shape and scale dependence of the leading 
order matter bispectrum. 

In reality, we can only measure the Fourier space density 


contrast at discrete grid points whose separation is deter¬ 
mined by the volume of the iV-body simulation, or the width 
of the survey window function. The bispectrum then can also 
be measured only for a finite number of triples (/ci,/c 2 ,fe). 
We can facilitate the comparison among different predictions 
for the bispectrum by flattening all possible triplets into a 
one-dimensional vector (x-axis), where we rank the triplet in 
row-major, ascending order with the condition ki > k 2 > ks. 
This allows us to look at all of the triangular configurations 
over a range of scales simultaneously, as well as compare 
bispectra on the y-axis. Each plot in the right column of 
Fig. 2 shows the bispectrum of all triangular configurations 
in a range of ki that includes the slice shown in the left 
panel. The x-axis shows the ranked (/ci,/c 2 ,fe) triplets, and 
the y-axis shows the amplitude of bispectrum. As one can 
see from the plots, a similar pattern block is duplicated for 
each value of the long-wavelength mode ki , as indicated by 
vertical dashed lines. We also show the configuration de¬ 
pendence by highlighting the limiting triangles in different 
colors. One can easily relate this plot to the contour plot in 
the left panel. For a given ki, the first-ranked triplet is at the 
bottom vertex (/c 2 — k^ — ki/2^ blue dots, folded triangle) of 
the triangular region in the left panels (k 3 /ki-k 2 /ki plane). 
The rank is followed by the next elongated triangle (green 
dots) then moves rightward to the isosceles triangles (red 
dots). This process continues until the rank reaches most- 
squeezed triplet then walks along the upper-most horizontal 
line to the equilateral triplet. Depending on the scales (ki 
values) that we are looking at, the configuration dependence 
of the bispectrum is different. 

In our analysis for this paper, we shall mainly use this 
flattened representation to study the configuration and scale 
dependence of bispectrum at all triangles. Although non- 
linearities change the exact shape of the bispectrum ‘chunk’ 
for a given ki , the overall structure remains unchanged. 


3 TRANSIENTS IN STANDARD 
PERTURBATION THEORY 

In this section, we review the standard perturbation the¬ 
ory (PT) formalism for calculating the tree-level bispectrum 
(Bernardeau et al. 2002). We then follow the approach of 
Scoccimarro (1998) and Crocce et al. (2006) to compute the 
correct growing and decaying modes for initial conditions 
set using the Zel’dovich Approximation and 2LPT initial 
conditions. 

One way of simplifying the single fluid equations is com¬ 
bining the density field 8{k^rf) and velocity divergence field 
0{k,r]) into a doublet as: 

^a{hri) = {S{k,ri),-0{k,ri)/n) (4) 

where the logarithm of the linear growth factor (D(r)) rel¬ 
ative to its initial value (Di): rj = \n{D{T)/Di) is the time 
variable. Here, 1-i = dlna/dr = aH is the reduced Hubble 
parameter that is defined with conformal time r = / dt/a{t). 
The equations of motion for Ta, combining the continuity 
equation, the Euler equation, and the Poisson equation, can 
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(b) Flattened bispectrum, with 5ki = 0.001 h/Mpc 
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(e) ki = 0.15 h/Mpc (f) flattened bispectrum, with 5ki = 0.005 h/Mpc 

Figure 2. Left: slices of the bispectrum at constant ki. The locations of different triangular conflgurations are highlighted. Right: the 
ffattened bispectrum with a given 5ki^ where the range of ki plotted includes the corresponding slice on the left. Vertical dashed lines 
indicate different slices of ki. Values of ki in several slices is given in units of h/Mpc. 




then be written as^: 

ri) + Q,ab'^b{k, 7 ?) 

= l‘':l{k,kuk 2 )^b{kuri)^c{k 2 ,ri)- ( 5 ) 

Here, 

Q.ab= _ 3^2 1/2 

^ As the initial conditions are set in the deeply matter-dominated 
epoch, and transient effects are most dominant at higher redshifts, 
our analysis throughout this paper is restricted to a ffat, matter- 
dominated (Einstein de-Sitter) universe. Later time dark energy 
will only marginally change the analysis as the PT kernels are 
insensitive to the background cosmology Bernardeau et al. (2002). 


is the linear mixing matrix, and 

fci, fe) = (2^)3fo(fc- fci - fc2)^ (7) 

^il\ik,ki,k2) = (27r)^fe(fe- fei - k2)’^^^^ (8) 

encode the non-linear interactions. Note that 
7121(^5^15^2) = 7^^2(^5 ^2, and the integration 

J dPki/{ 27 ^)^ over the repeated wavevectors {ki and 
fe) are implied. 

A Laplace transform of Equation ( 5 ) in the variable 77 
leads to: 

UJ^a{k, to) - (j)a{k) + VLah^h{k, Uj) ( 9 ) 

^c-|-ioo 7 

='yabc(^’^l’^2) / ^'Pb(kl,CUl)9c(k2,CU - CUl) , 

J c—ioo 


MNRAS 000, 1-?? (0000) 





















































































































Toward Accurate Modeling of the Nonlinear Bispectrum 5 


where 


(t)a{k) = ^a{k,r] = 0) 


( 10 ) 


is the density and velocity field at the initial time (i.e. the 
time at which we start the simulation). The linear part of 
Eq. (9) can be written in terms of a matrix aab’- 

^ab (^) ~ ^^ab T ^ab (^^) 


CTabicj) 


1 

(2u; + 3)(a;- 1) 


2,00 T 1 
3 


2 

2u ’ 


( 12 ) 


and, finally, we obtain the formal solution through an inverse 
Laplace transformation as: 


'^a{k,ri) = gab{ri)4)b{k) (13) 

+ / dr]'gab{ri-r]')'li%{k,ki,k2)^c{ki,r]')^d{k2,r]'). 

Jo 

Here, 

nc+ioo 1 

gabin) = / ^aabHe^^ (14) 

J c—too 27rz 



■3 

2 ■ 

e-3v/2 

■-2 2 

5 

3 

2 

5 

CO 

1 

CO 


is the Green’s function which describes the time evolution of 
the linear perturbations. The first term proportional to oc 
a is the growing mode and the second term proportional to 
g-3T7/2 ^ decaying mode. 

Eq. (13) is the formal solution of the fluid equations, 
and starting points of various of renormalized perturbation 
theories. Eor our purposes, however, it suffices to find the 
perturbative solution that follows. We perturbatively ex¬ 
pand the solutions at initial and final time (0a and Ta, re¬ 
spectively) in terms of the linear (Gaussian) density 8{){k) 


field at initial time as: 

0,(fc) = y34"^(fc) 

n=l 


(15) 

00 

'^aik,g) = y^^'i'*^(fe,r?) 

n — 1 


(16) 

where 

cPi^\k) = [dn]nTi'^\kl,- 

•• ,kn)6o{ki) ■ ■ ■ 6o{kn) 

(17) 


• • , kn,r])So{ki) • • • So{kn) 

(18) 


and [fcjn = dnik — ki — ■ — — kn). Here, are the n-th 
order kernels of the initial density and velocity fields, and 
are the n-th order resulting kernels at a later time 77 . 
Again, the integration f d^fc/(27r)^ over the repeated fe’s is 
implied. By using this ansatz, we calculate the kernels 
at the final redshift from the initial kernels Ta^^: 

= gab{ri)T^’^'> (19) 

+ E /" Sabig - 

m=l '^0 

This recursion relation is the starting point of the transient 
analysis. 


3.1 Growing mode solutions of Standard 
Perturbation Theory 

In the real universe, the decaying modes decay away at very 
early times, and the density field is at its fastest growing 


mode for any reasonable choice of the initial redshift. There¬ 
fore, the initial conditions for ideal iV-body simulations must 
be set at an early enough time such that PT is valid on 
the smallest resolved scales, and set by the fastest growing 
modes at all relevant orders in PT. 

In this section, we find the fastest growing mode so¬ 
lution perturbatively. We first find the linear growing mode 
(l>^a\k) — 5o(^)(l, 1) SO that the decaying part of the Green’s 
function Eq. (14) vanishes. Then, the second order kernels 
at time 77 are given by Eq. (19) using the linear grow¬ 
ing kernel of — e^(l, 1). While the time integral in 

Eq. (19) gives rise to the slowly growing mode (oc e^) and 
the decaying mode (oc those are exactly canceled 

by the initial fastest-growing kernels Ta‘^^ — {F 2 ^\g^ 2 ^) 
defined by the kernels in standard PT (Bernardeau et al. 
2002 ) and we are left with the fastest-growing solution of 
^ 2 ^^)- The same is true for all the higher 
order kernels, and the fastest growing solutions are 

■^f‘\k,'n) =i2nf5D{k-ki - K) 

X Tjrfki,.., kr,,g)5o{ki) ■ ■ ■ 5o{kn), (20) 


where 

(FlC\Gf^) ( 21 ) 

are the standard PT kernels. In particular, the second order 
kernel which appears in the leading-order matter bispectrum 
in Eq. (2) is given by 


J’P(fcl,fc 2 ) 


5 fci • fe / 1 1 \ 2 (fci • fc2)2 

2 Gi ^ 2 / 7 kfk^ 


( 22 ) 


and the velocity kernel is given by 


Gf\ki,k2) 


3 ki ■ k2 / 1 1 \ 4 (fei • fe)^ 

7”^ 2 Gr 7 fcjfci 


(23) 


When the single fluid approximation is valid (at early 
times and on large scales), ideal iV-body simulations must 
reproduce the density and velocity fields in Eq. (20) with the 
growing mode kernels Eq. (21). To satisfy this condition, we 
must start simulations with initial density and velocity fields 
satisfying Eq. (20). 


3.2 Transients from Initial conditions of iV-body 
simulations 

Generating the fastest-growing initial conditions at all or¬ 
ders in perturbation theory is quite non-trivial. This is in 
part because we simulate A’-body dynamics with the motion 
of matter partieles. That is, the initial conditions for A-body 
simulations are set by position and velocity of each particle 
(rather than set by density and velocity gradient field as is 
demanded by the PT analysis). Due to the nonlinear relation 
between particle’s position and the resulting density field, 
the initial conditions can only satisfy the fastest-growing 
condition perturbatively, in an order-by-order manner. As a 
result, when generating initial conditions with the growing 
mode at a given order in (5o, it is inevitable that spurious 
decaying modes will be excited at higher orders. 

In this section, we study the decaying modes for two 
of the most widely used initial condition generators: first 
order (Zel’dovich approximation) and second order (2LPT) 
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Lagrangian perturbation theory (LPT). For each initial con¬ 
dition generating scheme, we calculate the initial density and 
velocity kernels Ta^\ and find the evolved kernel at the fi¬ 
nal redshift by using Eq. (19) in order to study the resulting 
transients in the matter bispectrum as a function of initial 
and final redshifts as well as triangular configuration and 
scales. 

LPT describes the evolution of the matter density field 
by displacement ^ between each matter particle’s initial 
Lagrantian position q and final position x\ 

x{q,ri) = q+'^[q,ri), (24) 

then the peculiar velocity field is given by 

v{q, ri) = ^ = ^'{q, 7?). (25) 

Note that x is the comoving coordinate, and prime denotes 
the time derivative with respect to the conformal time r. 
When it comes to generating initial conditions of Wbody 
simulations, q refers to the position of particles on a regular 
grid, and x refers to particle’s initial position. 

For a given displacement field we find the Fourier 
space density contrast and velocity gradient in terms of the 
displacement field as 

S{k,r]) = J (f X6{x)e~^^'^ 

= E ^ [-ik ■ ^{q, vT , (26) 

n=l 

0{k,r]) = y d®gJ(9,7?)Vx •(27) 

where we calculate the Jacobian J(q,r]) by using mass con¬ 
servation J{q,r]) = \d^x/d^q\ = (l + (5(a?))“^ = \Sij + 

Here, comma stands for the derivative with respect to the 
Lagrangian coordinate. As the cosmic matter field on large 
scales is irrotational, we introduce the scalar potential such 
that Tj = In n-th order LPT, the scalar potential 0 

is expanded up to n-th order in the linear density contrast 
Jo (A;) (Bouchet et al. 1995). 


3.2.1 ZeVdovich approximation 

In the Zel’dovich approximation, the scalar potential 0 is 
given by 

(28) 

and from Eqs. (26)-(27) we find the linear order initial ker¬ 
nels is 

r,^^’(') = (i,i). (29) 


This means the Zel’dovich approximation indeed gives rise 
to the growing mode solution in linear order. On the other 
hand, the second order initial kernels are 


^ZA,(2) 

> a 


/(fc-fci)(fc-fe) k^{ki-k2)\ 

V 2fc2fc2 ’ 2klkl ) ’ 


(30) 


which are different from the fastest-growing mode kernels of 
PT in Eq. (21). Therefore, the transients in the Zeld’ovich 
approximation start to appear from second order, and affect 


the subsequent evolution of the density and velocity fields 
at second order as 

F^^’^^\ki,k2,v) = e'^^F^^\ki,k2,v) 

+ (iij{k^)Uj{k2) - I) (^e’' - (31) 

G 2 ^'^^\ki,k 2 ,r]) = e^'^G^ 2 ^\ki, fe,ry) 

+ (t,j{h)Uj{k2) - I) (^e’' + . (32) 


Here, tij{k) is the longitudinal operator which relates the 
density field and the tidal field Sij{k) 


,(fc) = 


fc 2 


-f S{k) = Uj{k)5{k), 


(33) 


which means that the Zel’dovich approximation contributes 
an extra tidal field to the iV-body simulation at second order. 

The extra tidal field then yields a transient in the lead¬ 
ing order matter bispectrum, and the slowest-decaying tran¬ 
sient is given by 

B^^{kl,k2,k3)-B{kl,k2,k3) 

3 

(mi 2 - 1) PL{ki)PL{k 2 ) + (2 cyclic), (34) 


3 ^ 
5 [Di 


with /j^ij = ki • kj being the angular cosine between two 
wave vectors. With Eq. (19), one can also show that, 
with Zel’dovich initial conditions, the slowest-decaying (or, 
longest-lasting) transients at all higher orders (n > 2) are 
suppressed only by with respect to the fastest-growing 
modes. 


3.2.2 2LPT 

In 2LPT, the scalar potential is given by 
0(fe) = ^ + {2nf5^ik-ki2 p\~,t'J^^ 5o{ki)Soik2), 


14fc2 


(35) 


with ki 2 = ^1 + ^ 2 , which reads, from Eqs. (26)-(27), 

7 ;^™) = ( 1 , 1 ) 

7;2LPT,(2) ^ 


(36) 

(37) 


Then, as shown in Sec. 3.1, the density and velocity fields in 
2LPT remain as the fastest-growing modes, and an iV-body 
simulation with 2LPT initial conditions correctly reproduces 
the leading order bispectrum. 

On the other hand, 2LPT initial conditions differ from 
the fastest-growing modes at third order. This generates spu¬ 
rious slowly growing modes and decaying modes. Again, us¬ 
ing Eq. (19), one can show that, for 2LPT initial conditions, 
the slowest-decaying modes are suppressed by two powers in 
the linear growth factor (e“^^) with respect to the fastest 
growing mode at n(> 3)-th order. Therefore, we expect the 
following time dependence of transient: 


A:2, ks) - H(/ci, /c2, fe) oc 



(38) 


and to contribute on smaller scales than the Zel’dovich case 
because the transient only comes from next-to-leading order 
and beyond. Relative to the leading-order, fastest growing 
mode, the transient decays faster than the Zel’dovich case 
as (D/Di)-^. 
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4 TRANSIENTS IN SIMULATIONS 

In this section, we study the error induced by transients 
numerically with a suite of simulations with different ini¬ 
tial redshifts, using both the Zel’dovich approximation (ZA) 
and 2LPT for setting the initial positions and velocities of 
particles. 

We examine a set of twelve simulations with different 
initial redshifts, but with the same random seed to set ini¬ 
tial conditions. This ensures that within each set, any differ¬ 
ences we measure between the simulations are due solely to 
the initial conditions generator. The initial conditions in six 
of the runs were generated using ZA, and in the other six 
we used 2LPT. We use the following initial redshifts in each 
case: zinit = 400, 300, 200,100, 50, 25. We use the 2LPT sim¬ 
ulation started at zinit = 400 as our reference simulation to 
which all other runs are compared. That is, this simulation 
defines our fiducial statistics, and 

We run the simulations using the TreePM part of the 
publicly available version of Gadget 2 code. Each simula¬ 
tion is run with 512^ particles in a 200 Mpc/h box and 
particle-mesh Fourier grid iVgrid = 1024^. The background 
universe is ACDM with the following cosmological parame¬ 
ters (Qa = 0.73, Qrn = 0.27, h = 0.7), and the initial density 
field is calculated from the primordial curvature power spec¬ 
trum with spectral index Ug = 0.96, and normalized so that 
the r.m.s. fluctuation smoothed by a spherical top-hat filter 
with radius 8 Mpc/h is crs = 0.8 at present (z = 0). To min¬ 
imize the effect from accumulated numerical noise, we have 
required stringent accuracy of N-body simulation by setting 
parameters of Gadget 2 following the choice of Crocce et al. 
(2006): ErrTolIntAccuracy = 0.01, MaxRMSDisplacement- 
Fac = 0.1, MaxSizeTimestep = 0.01, ErrTolTheta = 0.2, 
ErrTolForceAcc = 0.002. We calculate the cloud-in-cell den¬ 
sity from the output particle distribution on a 1024^ grid at 
z = 6, 5, 4, 3, 2 and 1 and compute the power spectrum and 
bispectrum at each redshift. 


4.1 Transients in the matter power spectrum 

Figure 3 shows the fractional error in the measured power 
spectrum, as compared to the fiducial power spectrum, at 
z = 1,2,3, and 6, from 2LPT (solid) and ZA (dashed) ini¬ 
tial conditions. The various initial redshifts are shown in 
different colors. As expected, the transients from ZA initial 
conditions are in general larger than 2LPT on all scales and 
redshifts, apart from on very small scales at z = 6, and the 
error is larger for smaller initial redshifts in all cases. While 
the 2LPT power spectrum transients decay quickly, the ZA 
transients remain significant, even on relatively large scales, 
until z = 1. These findings are consistent with PT predic¬ 
tions for transients in the matter power spectrum, and also 
with previous studies Crocce et al. (2006); Jeong (2010). 


4.2 Transients in the matter bispectrum 

Figure 4 shows the fractional error in the measured bispec¬ 
trum compared to the fiducial bispectrum, again a,t z = 
1, 2, 3, and 6, from 2LPT (top) and ZA (bottom) initial con¬ 
ditions. Note the different y-ranges in the two panels. As 
we found with the power spectrum transients, the error is 
larger for smaller initial redshifts. The bispectrum transient 


error from the ZA simulations is significantly larger on all 
scales than in the 2LPT simulations, because the ZA tran¬ 
sients appear at tree-level as opposed to 1-loop in PT. We 
also note that the transient error is larger at higher red- 
shift in all cases, as expected from theory. We also show the 
transient effect in the corresponding reduced bispectrum 


Qi^ki , /C2 , ^3 ) 


R(/Cl,/C2,fe) 

P{ki)P{k2) + P{k2)P{k3) + P{k3)P{ki) 

(39) 


in Figure 5. As expected, when removing the effects due to 
the transients in the power spectrum, the transient effects in 
the reduced bispectrum are smaller, in particular for lower 
redshifts and smaller scales, than the transients in the bis¬ 
pectrum. The transients in 2LPT simulations, however, are 
much smaller than those in ZA simulations. 

In order to accurately model the nonlinear bispectrum 
numerically, we suggest using 2LPT initial conditions with 
an initial redshift zinit > 100. For Zi = 100, the error from 
transients is less than 1% for z < 3 and around 2% at z = 6 
for this box size and resolution. For the corresponding ZA 
simulation, the errors are around 8% for z < 3 and 10% for 
z = 6. Although the degree of the transient effect depends 
on other factors such as the resolutions of simulation, we find 
this rule of thumb also works for simulations with Lbox = 
1 Gpc//Z- and -/Vparticle — 512 . 

Another consideration is the shape-dependence of the 
bispectrum transients. Figure 6 shows the transient signal 
{AB — B — top) and induced error [AB/B^^] bot¬ 
tom) at ki = 0.7 h/Mpc (normalized to —1). On the left 
is the prediction from PT for the transients from ZA initial 
conditions. The middle column shows the signal and error 
measured from the the Zi = 200 ZA simulation at z = 1. 
The right column shows the same for the 2LPT simulation 
with Zi = 200. 

While Figure 4 clearly shows that the 2LPT transients 
have a much smaller amplitude than the ZA transients. Fig¬ 
ure 6 shows that they have very similar shape-dependence to 
ZA transients on these scales. The measured shape depen¬ 
dences of the transient signal (top panels) and error (bottom 
panels) agree reasonably well with the SPT prediction for 
both ZA and 2LPT ICs. The error induced from transients 
is greatest in the equilateral configuration, and smallest in 
the elongated configuration for both ZA and 2LPT initial 
conditions. We find that the shape dependence of the tran¬ 
sients is independent of redshift in both ZA and 2LPT cases; 
it is only the amplitude of transients that is reduced at lower 
redhsifts because the decaying modes fade out as the simu¬ 
lation proceeds. 


5 MODELING NONLINEARITY IN THE 
MATTER BISPECTRUM 

We now test the accuracy of various theoretical templates for 
the nonlinear bispectrum against the measured matter bis¬ 
pectrum from A-body simulations. For comparison, we use 
the result of zinit = 400 (2LPT) simulation, as the nonlin¬ 
ear matter bispectra measured from other 2LPT simulations 
converge to it. 

On large enough scales, the leading order (tree level) 
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(b) Transient effect at 2 ; = 2 



<1 



(d) Transient effect at 2 ; = 6 


Figure 3. Error induced by transients in the measured power spectrum at redshifts 2 ; = 1,2, 3, 6 using ZePdovich and 2LPT initial 
conditions with different initial redshifts. In each plot, the left column shows the results from the large-box simulation and the right 
column shows the small-box simulation. The top panels of the plots show the error in simulations using 2LPT initial conditions and the 
bottom panels show the error from ZePdovich initial conditions. 


prediction in Eq. (2) is enough to model the matter bis¬ 
pectrum. Modeling the matter bispectrum on smaller scales 
where nonlinearities become stronger, however, requires 
more elaborate calculation. On quasi-linear scales, we ex¬ 
pect the next-to-leading order (one-loop) bispectrum in PT 
to model the nonlinear evolution reasonably well. The next- 
to-leading order matter bispectrum is given by four addi¬ 
tional one-loop terms (Bernardeau et al. 2002): 

fc2, fca) fc2, fca) + {kiM, ka) 

+B^^^^'^\kuk2,k3) + B^^^^’’\kuk2,k3) 
+B(“^)(fci,fc2,fc3), (40) 

where is the leading order bispectrum in Eq. (2) and 

B^^^^\ki,k2,k3) = 8 J - q) 

X k 2 + q)F^"\q - fei, -fe 2 - q) 

^PL{q)PL{\k,-q\)PL{\k 2 + q\) (41) 

B(“^Hfci,fc2,fc3) = 12Pz,(fc2)Bi,(fc3) f -^PLiq) 

X -q, -fe, -ks) + (2 eye.) (42) 

B('2"“)(fci,fc2,fc3) = 6PL(fci) J 0^PLiq)PLi\k2-q\) 

X Fppg., h 2 - q)F^^\-q,q-h 2 , -fei) + (5 eye.) (43) 


B<'"®'’Pfcl,fc2,fc3) = 6PL{kl)PL{k3)F^^\kl,k3) 

are non-linear corrections with being symmetric kernels 
encoding the fastest growing mode of the nonlinear density 
evolution (see, e.g. Sec. 3). Note that the expression above 
is completely determined by the given linear matter power 
spectrum PL(k). 

On smaller scales where PT breaks down due to strong 
nonlinearities, various phenomenological models have been 
proposed as alternative ways of predicting the nonlinear be¬ 
havior of the matter bispectrum. Here, we consider three 
such models (PT+, GM model Gil-Marm et al. (2012), and 
SG model Scoccimarro & Gouchman (2001), see below for 
details) for the nonlinear matter bispectrum. 

The simplest model that we consider is PT+, where 
we attribute the entirety of the nonlinear evolution of the 
matter bispectrum to the nonlinear matter power spectrum. 
This leads to an expression that is identical to Eq. (2), but 
uses the measured (fully-nonlinear) power spectrum from 
the iV-body simulation in place of the linear power spectrum: 

B(fci,fc2,fc3) = 2Fppfci,fe2)P(fci)P(fc2) + (2 eye.) (45) 

This model is motivated from the expression of the 
one-loop bispectrum in Eq. (40), which can be written 
as P“°°P(fcl,fc2,fc3) = 2P2^®^(fci,fc2)Plloop(fcl)Plloop(fc2) + 
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(a) Transient effect at 2 ; = 1 (b) Transient effect at 2 ; = 2 




(d) Transient effect at 2 : = 6 


Figure 4. Error induced by transients in the measured bispectrum at redshifts 2 : = 1, 2, 3, 6 using ZePdovich and 2LPT initial conditions 
with different initial redshifts. The top panels of the plots show the error in simulations using 2LPT initial conditions at various initial 
redshifts and the bottom panels show the error from ZePdovich initial conditions. The different line colors indicate different initial 
redshifts for the simulations. Vertical dashed lines indicate slices of fci, and the ki values of several slices are indicated along the a^-axis 
in units of h/Mpc. 


.z =2 
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(a) Transient effect at 2 ; = 1 


(b) Transient effect at 2 ; = 2 







fci = 0.31 


fcj = 0.75 



(c) Transient effect at 2 ; = 3 (d) Transient effect at 2 ; = 6 

Figure 5. Same as Figure 4, but for the reduced bispectrum, Q. 
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0.0 

- 0.1 

- 0.2 

-0.3 

-0.4 

-0.5 

- 0.6 

-0.7 

- 0.8 

-0.9 

- 1.0 


Figure 6. Shape dependence of the transient signal and error in the measured bispectrum at ki =0.7 h/Mpc at 2 : = 1 from Zel’dovich 
and 2LPT initial conditions compared to the SPT prediction. The top panels show the transient signal, which is highest in the squeezed 
and elongated configurations. The bottom panels show the percent error induced by transients, which is highest in the equilateral 
configuration. Each slice is normalized by its minimum value, so the range is from —1 to 0 in each case. 


2(eye.)+(irreducible terms), with the one-loop power spec¬ 
trum Piioop(^)- In the language of renormalized perturba¬ 
tion theory Crocce & Scoccimarro (2006), PT+ corresponds 
to fixing the ‘vertex’ (the kernel while using the renor¬ 
malized propagator (power spectrum P{k)). 

PT+ is a simplified version of the fitting formula pro¬ 
posed by Scoccimarro & Couchman (2001) (hereafter, SC), 
where the ^ 2 ^^ kernel in Eq. (45) is replaced by the effective 
kernel as follows. 


kernel are replaced by tilded functions as 


d{n, k) = 


(50) 

1 + (gai)’*+“2 

b{n, k) = 

1 + 0.2a3(n + 3)g’*+®(ga7)"+®+“® 

(51) 

1 + (^a7)^+^-^+“8 

c(n, k) = 

1 + 4.5a4/[1.5 + {n + 3)^](ga5)'*+^+“® 

(52) 

1 + (^a5)’^+^-^+'^9 


F2{ki,kj) = -a(ni,ki)a{nj,kj) 

+ icos(6>ij) + h) b(rn,ki)b(nj,kj) 

2 2 

+ -COS {dij)c(ni,ki)c(nj,kj). (46) 

Here, Oij is the angle between two wavevectors fe and kj, 
and reduces to the PT kernel in Eq. (22) when a = b = 
c = 1. The three functions a(n, /c), 5(n, /c), c(n, k) are given 
as functions of the wavenumber k and the effective slope of 
the linear power spectrum n = d\nPL{k)/d\nk as 


a{n, k) 
6(n, k) 
c{n, k) 


1 + a^^z)[0.7Q3{n)Y^^{qair+^^ 

1 + 

l + 0.2a3(n + 3)g^+^ 

1 + qn+s.5 

1 + 4.5ci4/[1.5 + (71 + 3)"^] (^0,5)’^”*"^ 
1 + 


(47) 

(48) 

(49) 


where q = k/kn\ with nonlinear scale kn\ defined by the 
wavenumber at which the dimensionless power spectrum 
A‘^{k) = k^PL{k)/{27r‘^) becomes unity A^(A:ni) = 1, and 
varies with redshift as A^ni oc (1 + z). Q 3 (n) = (4 — 2’^)/(l + 
2^^^^). SC fit the parameters ai through ae using the bis¬ 
pectrum measured from P^M simulations Jenkins et al. 
(1998) with 256^ particles in a 240 Mpe/h box to find that 
ai = 0.25, a 2 = 3.5, as = 2, a4 = 1, us = 2, ae = —0.2. 

Recently, Gil-Marm et al. (2012) (hereafter, CM) pro¬ 
posed extending the SC model further with 3 additional pa¬ 
rameters, where a(n,/c), b{n,k) and c(n,/c) in the effective 


The definitions for n, /cni, and Q 3 are the same as above, 
except in this model n{k) is smoothed so the oscillatory fea¬ 
tures from the Baryon Acoustic Oscillations (BAO) are re¬ 
moved (further discussion of this can be found in Gil-Marm 
et al. (2012)). Then, CM fit the parameters ai to ag us¬ 
ing the measured matter bispectrum from two sets of iV- 
body simulations, one with L = 2.4 Gpc/h (768^ parti¬ 
cles), and one with Lb = 1.875 Gpc/h (1024^ particles), 
to find the best-fit parameters of ai = 0.484, a 2 = 3.740, 
as = —0.849, a4 = 0.392, as = 1.013, ae = —0.575, 
ar = 0.128, ag = -0.722, ag = -0.926. 

Note that the conditions on which the fitting formulae 
are developed are quite different. While SC have used the 
matter bispectrum at all triangular configurations (a total of 
a million triangles for /c < 2.3 h/Mpe) at redshift z = 0 and 
z = 1, CM have only used certain configurations (^ 12/77 = 
0.1, 0.2, •••, 0.9 and = 1.0, 1.5, 2.0, 2.5, for fe < 

0.4 Mpe/h) at four different redshifts z = 0, 0.5, 1, 1.5. 
Moreover, they are subject to different levels of transient 
errors as the ACDM simulations used by SC are generated 
at Zinit = 30 with the Zel’dovich approximation, and the 
simulations used by CM are generated by 2LPT at zinit = 19 
(2.4 Gpc/h box) and zinit = 49 (1.875 Gpc/h box). 

To test the validity of each of the phenomenological 
formulae outside of the dynamical (Eourier) ranges and red¬ 
shift ranges at which the models are fitted for, we com¬ 
pare the model predictions to measured bispectra at red- 
shifts 1 < z < 6 in both large-scale (Lbox = 1 Gpc/h) and 
small-scale (Lbox = 200 Mpe/h) simulations. We also com¬ 
pare these to the predictions from tree-level PT, PT+, and 
1-loop PT. Eor both the SC and CM models, we calculate 
the effective slope of the linear power spectrum n{k) from 
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k, [h/Mpc] 


Figure 7. The bispectrum at 2 ; = 0 measured from the 1 Gpc/h 
simulation, compared to the SC (green), CM (blue), PT (grey), 
and PT+ (red) models. The top plot shows the full flattened 
bispectrum, and the bottom plot shows only the squeezed and 
equilateral triangles as a function of ki. 

the BAO-smoothed linear power spectrum which provides 
better agreement than the n{k) including BAO. 

First, we check that we can reproduce the results from 
Gil-Marm et al. (2012) at 2 : = 0 on large scales. We expect 
that the GM model is the most accurate here, because their 
model parameters were fit from a similar simulation. Fig¬ 
ure 7 shows the measured bispectrum at 2 ; = 0 from the 1 
Gpc//i simulation, compared to the SG and GM models, as 
well as PT and PT+. The top plot shows the full flattened 
bispectrum, and the bottom plot shows only the squeezed 
and equilateral triangles as a function of ki . We can see that 
the GM model indeed gives an improved agreement to the 
measured bispectrum over the SG model on these scales at 
2 : = 0. 

Next we look at 2 ; = 1 in both the 1 Gpc//i and 200 
Mpc/h simulations. Figure 8 compares the measured bis¬ 
pectrum to the models at 2 ; = 1. The left column shows the 
results from the large-box simulation, and the right shows 
the small-box simulation. Note that the ki range in the bot¬ 
tom plots overlaps, but we plot the results from the two 
simulations separately for clarity. Also note that the mea¬ 
sured signal for squeezed triangles from the two simulations 
is slightly different due to the different resolutions. From 
these plots, we see that the GM model still gives the best 
fit to the measured bispectrum on the scales shown in all 
triangular configurations, except for squeezed triangles on 
very small scales, where the SG model fits slightly better. 

However, we find at higher redshifts, the GM model is 
less accurate than the SG model. Figs. 9, 10, and 11 show the 
measured bispectrum at 2 ; = 2, 3, and 6, respectively, from 
both simulations compared to the models. In all of these 
cases, the SG model gives the most accurate prediction on 
all scales and in all triangular configurations. The GM model 
over-predicts the bispectrum slightly in all cases. On large 
scales (left panels), the predictions from perturbation theory 
(1-loop and PT+) give reasonably accurate predictions. 

It is clear from these comparisons that one must care¬ 


fully consider the range of validity of a given model when 
predicting the nonlinear bispectrum. At 2 ; < 1, the GM 
model undoubtedly gives the best prediction, but at higher 
redshifts, it is better to use the SG model. At high redshift 
( 2 ; > 3), the PT models are sufhcient for modeling the non¬ 
linear bispectrum. Further work is needed to develop a more 
reliable model of the nonlinear bispectrum on all scales and 
redshifts. 


6 CONCLUSION 

As the non-linear growth of structure generally causes 
the cosmic density field to become non-Gaussian, study¬ 
ing higher-order statistics beyond the two-point correlation 
function will open up a new avenue to exploit a wealth of 
cosmological information. In particular, with a large-volume 
coverage and high number density of galaxies, current and 
future galaxy surveys promise ever-more accurate measure¬ 
ments of the galaxy three-point function and bispectrum. 
These measurements of three-point correlation functions can 
give us a better understanding of the growth of structure, 
galaxy bias, and primordial non-Gaussianity. 

Extracting information from the galaxy bispectrum re¬ 
quires accurate modeling of the bispectrum of the dark- 
matter density field. We have shown that when studying 
the matter bispectrum numerically, it is important to un¬ 
derstand the effects of transients from initial conditions in 
simulations. When generating the initial conditions using 
Zel’dovich approximation, the transient of the matter bis¬ 
pectrum shows up at leading order in perturbation theory 
and affects all scales. In contrast, the transient of the mat¬ 
ter power spectrum shows up only at higher-order and af¬ 
fects mainly small scales. In order to simulate nonlinearity 
in the matter bispectrum correctly, a 2LPT (second-order 
Lagrangian perturbation theory) initial conditions genera¬ 
tor must be used. To further reduce the effects of higher- 
order transients, we recommend using 2LPT initial condi¬ 
tions with Zinit > 100. This requirement is based on our 
suite of simulations, in order to maintain a sub-percent level 
accuracy of the matter bispectrum on scales k < 1 h/Mpc 
at z < 3. 

Another challenge is theoretical modeling of the nonlin¬ 
ear behavior of the bispectrum beyond tree-level perturba¬ 
tion theory. We have shown that analytical calculations of 
the next-to-leading order, one-loop bispectrum models the 
non-linearities in the matter bispectrum quite well on quasi- 
nonlinear scales, in particular at high redshifts. The fact 
that this analytical calculation of the one-loop bispectrum, 
which does not contain any free parameters, can model the 
non-linearities in the matter bispectrum on quasi-nonlinear 
scales is very encouraging for theoretical modeling. More sys¬ 
tematic and in-depth studies, parallel to Jeong & Komatsu 
(2006) but for the matter bispectrum, are needed in order 
to assess the accuracy of perturbation theory modeling at 
high redshifts in the quasi-nonlinear regime. 

On even smaller, non-linear scales, we found that the 
fitting formula from Gil-Marm et al. (2012) is reliable at 
low redshift {z < 1), but at higher redshifts, Scoccimarro & 
Gouchman (2001) gives the best agreement with the bispec¬ 
trum measured from simulations. We caution the readers to 
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Figure 8. The bispectrum at 2 ; = 1 measured from the large-box (left) and small-box (right) simulations, compared to the SC (green), 
GM (blue), PT (grey), and PT+ (red) models. The top plot shows the full flattened bispectrum, and the bottom plot shows only the 
squeezed and equilateral triangles as a function of ki. 




Figure 9. Same as Fig. 8 but for z = 2. 


critically check the range of validity of these formulae before 
applying them. 
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